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INTRODUCTION 


Background 

In  estimating  guidance  error  dispersion  of  laser  guided  projectiles, 
one  quickly  comes  to  recognize  the  importance  of  laser  spot  motion  as 
a primary  external  error  source.  Over  a period  of  several  years  the 
author  and  his  colleagues  have  been  involved  in  making  such  estimates 
for  the  Copperhead  (CLGP)  projectile  and  its  predecessors.  In  this 
connection  I have  frequently  emphasized  the  importance  of  this  error 
source  and  its  significance  for  the  guidance  accuracy  of  several  flight 
vehicles.  References  [1]*  and  [2]  illustrate  this  emphasis.  As  a con- 
sequence of  the  importance  attached  to  laser  spot  motion  and  related 
aspects  of  the  laser  signature  such  as  pulse  dropout  and  laser  spillover, 
considerable  efforts  have  been  spent  studying  tracking  error  for  various 
designators  and  trackers.  Some  analytic  models  of  the  tracking  error 
are  given  in  References  [3]  and  [4].  In  particular.  Reference  [4]  describes 
the  statistical  techniques  used  to  produce  an  analytic  model  descriptive 
o>  the  tracking  error  of  the  Ground  Laser  Locator  Designator  (GLLD) 
observed  during  the  CLGP  OT-1  Tracking  Tests.  Each  of  the  above  models 
are  descriptive  of  the  performance  of  a particular  combination  of  tracking 
device  and  tracker  operator  faced  with  a particular  tracking  problem. 


Square-bracketed  numbers  refer  to  cited  literature. 


[1] 


[2] 


Memorandum  for  Record,  AMSAR-SAM,  23  July  1975,  subject: 
Guided  Projectile  Effectiveness  Study. 


Army-Navy 


Schlenker,  G.  J.  and  Heider,  R.  D.,  Distribution  of  Angle  of  Obliquity 
of_  Laser-Guided  Projectiles  With  Respect  to  the  Target  at  Imparl-. 

DRSAR/ SA/N-51 , (AD  A032683) , US  Army  Armament  Command,  Rock  Island,  IL, 
August  1976.  ’ ’ 


[3]  Pastrick,  H.  L.  and  Hollman,  H.  C.,  Analysis  and  Digital  Simulation 
jlp/telg  f<?r  C^GP: — Martin  Marietta  Aerospace  Design.  Report  No.  RG-75-29, 
(.Appendix  G,  Pulse  Dropout  and  Pulse  Dither  Subroutines),  US  Army 
Missile  Command,  Huntsville,  AL,  Dec  1974. 

[4]  AMSAR-SA,  10  Dec  74,  subject:  Reduction  and  Analysis  of 

CLGP-OT-1  Tracking  Data. 


11 


Particular  models  of  this  type  are  statistical  summaries  suitable  for 
a digital  computer  simulation  of  the  process  of  spot  motion  as  seen  by 
the  designator,  which  can  be  used  as  input  to  terminal  guidance  simulations 
and  laser  target  reflectance  models.* 

An  accurate  characterisation  of  laser  spot  motion,  whatever  its  source, 
is  also  important  because  (stochastic)  process  dynamics  are  also  pertinent 
to  the  problem  of  laser  spillover  and  its  consequence  for  guidance  accuracy. 
In  Reference  [5],  which  addresses  this  problem,  spot  motion  was  treated 
as  a second-order  stochastic  process. 

In  an  effort  to  standardize  laser  designator  testing  procedures, 
increased  attention  (Reference  [6])  has  been  given  to  various  phenomena 
which  give  rise  to  laser  spot  motion:  device  jitter,  atmospheric  propagation 

effects  (scintillation  and  beam  steering),  and  designator  tracking  error. 

In  a field-test  environment,  these  phenomena  all  contribute  and  cannot 
be  isolated  by  direct  measurement.  Consequently,  some  of  the  past  analyses, 
which  have  considered  spot  motion1 * * * 5 6  and  "tracking  error"  as  interchangeable 
terms,  have  been  somewhat  indiscriminate.  Pragmatically  this  lack  of 
attention  to  the  components  of  spot  motion  may  not  have  been  of  great 
significance  in  that  the  joint  process  is  what  is  of  concern  in  guided 
projectile  accuracy  analysis.  However,  a closer  attention  to  the 
components  of  spot  motion  does  produce  clarity  in  understanding  the 
disparate  phenomena. 


* In  studies  done  at  ARMCOM  the  laser  signatures  seen  by  the  projectile 
seeker  are  generated  by  a laser  target  reflectance  model  using  a T55 
tank  as  the  target.  Validation  of  this  model  is  contained  in  the 
following  report:  Beard,  J.,  Rice,  D.,  and  Ladd,  D.,  Target  Reflection 

Illumination  Model  With  Second  Order  (TRIMS') , (CONF)  , ERIM  Report  No. 
192200-2-F,  Environmental  Research  Institute  of  Michigan,  Aug  1975. 

Updated  documentation  of  this  type  of  model  will  be  included  in  a report 
in  preparation:  Ladd,  D*,  A Multifaceted  Target  Signature  Model  With 

Digital  Imaging.  ~ 

[5]  Schlenker,  G.  J.,  "Proportion  of  Energy  Spilled  Over  a Target  During 
Tracking  With  a Laser  Designator  and  Implications  for  Terminal 
Guidance,  Systems  Analysis  Directorate  Activities  Summary  November 
-2=-9_7-6 * RePort  No.  DRSAR/SA/N-60 , US  Army  Armament  Command,  Rock  Island, 
IL,  December  1976. 

[6]  Minutes  of  a Meeting  at  White  Sands  Missile  Range,  NM,  STEWS-TE-AG, 

18  June  1976,  subject:  Uniform  Standards  for  Laser  Designator  Develop- 

ments. 
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Overview  of  the  Report 

This  report  reexamines  some  of  the  spot  motion  data  from  the  CLGP 
0T“1  Tracking  Tests  with  the  intention  of  identifying  and  statistically 
isolating  the  stochastic  component  attributable  to  human  tracking  error 
from  those  produced  by  other  mechanisms.  For  notational  convenience 
thruout  the  balance  of  this  report,  device  jitter  and  beam-steering 
effects  due  to  the  atmosphere  will  be  subsummed  by  the  term  "scintillation 
component."  This  simplification  is  considered  justifiable  because  of 
the  dominance  of  atmospheric  effects  in  the  example  treated  here.  In 
decomposing  laser  spot  motion  into  tracking  (t)  and  scintillation  (s) 
components,  one  exploits  the  fact  that  these  components  have  different 
dynamics  which  are  nearly  separable  in  the  frequency  domain. 

Much  of  this  report  is  devoted  to  deriving  mathematical  expressions 
useful  in  characterizing  the  stochastic  spot  motion  process  and  in  pro- 
ducing digital  implementations  for  simulation  purposes.  Mathematical 
results  are  generally  presented  in  maximum  generality  with  other  applications 
in  mind.  Numerical  examples  are  provided  thruout  the  report  to  give 
the  reader  a (hopefully)  better  sense  of  magnitude.  For  convenience  of 
exposition,  the  report  is  divided  into  the  following  six  chapters! 

Chapter  1.  Decomposition  of  Laser  Spot  Motion  into  Tracking  and 
Scintillation  components. 

Chapter  2.  The  Analytical  Autocovariance  Function  for  Laser  Spot  Motion 
Developed  From  the  Spectral  Density. 

Chapter  3.  Digital  Computer  Implementation  for  the  First-Order  Component 
of  the  Spot  Motion  Error. 

Chapter  4.  Digital  Computer  Implementation  for  the  Second-Order  Component 
of  the  Spot  Motion  Error. 

Chapter  5.  Autospectra  of  the  Digital  Error  Processes. 

Chapter  6.  Spectral  Moments  and  Parameter  Estimation. 

Computer  source  program  listings  for  pertinent  methods  are  presented 
in  the  Appendicies. 

Summary 

In  Chapter  1 it  is  argued  that  the  autospectrum  of  spot  motion  consists 
mainly  of  separable,  additive  components:  one  due  to  the  human  tracking 
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process,  possessing  approximately  second-order  dynamics,  and  the  other 

due  to  processes  having  approximately  first-order  dynamics  with  a significantly 

higher  crossover  frequency  than  the  human  tracker  component.  Expressions 

for  the  gain  constant  in  each  component  of  the  autospectrum  are  developed 

in  terms  of  the  variance  of  each  component.  A numerical  example  of  the 

analytic  autospectrum  is  compared  with  an  experimental  estimate  derived 

from  the  azimuthal  error  of  Run  75A  of  the  CLGP  OT-1  Tracking  Tests. 

Parameter  values  of  the  analytic  model  were  selected  by  trial  and  error, 
under  constraints,  to  produce  an  acceptable  Chebychev  fit  to  the  experimental 
autospectrum  in  the  sense  of  minimizing  the  maximum  residual. 

It  is  sometimes  useful  to  make  comparisons  between  analytic  models 
and  experimental  estimates  in  the  time  domain  as  well  as  in  the  frequency 
domain.  Altho  time-  and  frequency-domain  descriptions  are  equivalent, 
autocorrelations  are  often  more  intuitively  appealing.  For  example, 
the  time  interval  which  must  separate  two  time  series  segments  before 
they  are  effectively  uncorrelated  can  be  observed  immediately.  Chapter  2 
develops  an  expression  for  the  autocovariance  (and  autocorrelation) 
function  for  the  spot  motion  process  and  extends  the  example  of  Chapter  1 
by  comparing  the  autocorrelation  function  of  the  analytic  model  with  an 
experimental  estimate. 

Chapter  3 treats  the  digital  implementation  of  the  first-order  component 
of  spot  motion.  Implementation  is  achieved  by  low-pass  filtering  of 
gaussian  white  noise.  Expressions  are  developed  for  the  coefficients 
of  the  digital  filter  and  for  the  variance  of  the  input  noise  required 
to  produce  a given  output  process  variance. 

Chapter  4 parallels  the  developments  in  Chapter  3 while  treating 
the  second-order  component  of  the  spot  motion  process . This  development 
is  the  basis  of  the  principal  algorithm*  which  the  ARMCOM  Systems  Analysis 
Directorate  and  others  have  used  to  describe  laser  spot  motion  during 
the  past  four  years.  An  updated  version  of  the  computer  program  which 
we  use  to  describe  spot  motion  is  found  in  Appendix  A. 


* A special  case  of  this  algorithm  is  presented  in  subroutine  DITH  of 
Reference  [3], 
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Chapter  5 develops  the  autospectra  of  the  scintillation  (first-order) 
and  tracking  (second-order)  dynamic  components  as  implemented  digitally. 

The  autospectrum  of  each  component  is  separately  compared  with  its  con- 
tinuous-time (analog)  counterpart.  Some  sources  of  distortion  in  the 
digital  implementation  and  in  the  process  of  making  statistical  estimates 
of  autospectra  are  identified  and  are  quantified  using  the  first-order 
process  as  an  example.  These  results  suggest  the  means  to  minimize  spectral 
distortions . 

Chapter  6 derives  expressions  for  the  spectral  moments  of  each  of 
the  stochastic  components  of  the  spot  motion  process.  Additionally, 
statistical  estimators  are  derived  for  the  parameters  of  pure  first-order 
and  Butterworth  second-order  processes  which  exploit  the  spectral  moments. 
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CHAPTER  1 


DECOMPOSITION  OF  LASER  SPOT  MOTION 
INTO  TRACKING  AND  SCINTILLATION  COMPONENTS 


Let  the  variance  of  the  spot  motion  at  the  target  due  to  only  tracker 
2 2 

error  be  a in  (m  ) . Experience  with  the  TOW  and  GLLD  trackers  as  well 
as  airborne  man-operated  trackers  shows  that  tracking  error  exhibits 
approximately  second-order  dynamics. 

Let  the  tracking  process  be  second-order  Butterworth  so  that  the 
spectral  density  is  given  by 


H (V)  = 


1 + (v/vt) 


, 0 < V <°° 


with  gain  constant  A and  natural  corner  frequency  v . 
The  variance  or  total  power  of  the  process  is 


2 

a = 

t 


(v)  dv 


= A V 


df 

o i + tL 


(1.1) 


(1.2) 


(1.3) 


2 

a 


A V 7T  2 / 4 
t 


Whence 


2 

A = 4 cr 

JT  77  V 

A = 2 /T  a2 
t 


(1.4) 


(1.5) 


(m2/hz) 
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Normally,  the  corner  frequency  v is  determined  by  the  dynamics  of 

the  human  motor  system  and  lies  below  one  hertz. 

£ 

Scintillation  data  suggest  that  the  majority  of  scintillation  power 

lies  above  one  hertz  at  frequencies  to  which  the  human  cannot  respond. 

Consequently,  tracking  error  noise  and  scintillation  noise  are  substantially 

uncorrelated.  Data  on  intensity  of  scintillation  above  10  hertz  suggest 

this  process  has  nearly  first-order  dynamics. 

We  suppose  that  a first-order  process,  uncorrelated  with  the  tracking 

2 

process,  with  variance  O is  superposed  thereon.  The  spectrum  for  this 

s 

scintillation  process  is  given  by 


Hs  (V) 


B 

1 + (v/v  )2 

s 


with  corner  frequency  V . 

s 

Therefore, 


2 

a 

s 


H2  (v)  dv 

Jo  S 


a 


2 


s 


BV 

s 


df 

Jo  1 + f2 


(1.6) 


(1.7) 


(1.8) 


a2=  Bv  1 

s s 2 


Whence, 


s 


(1.9) 


(m2/hz) . 


(1.10) 


*Some  Power  spectra  for  percent  modulation  are  given  on  pp  224  of  [7],  Wolfe 

°f  Military  !nfrared  TschnnWv  Office  of  H„.l  R^d., 
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If  the  processes  are  superposed,  the  spectrum  of  the  joint  process 

is  just  the  sum  of  the  component  spectra.  Letting  the  joint  autospectrum 

2 * 
be  H (V), 

H2(V)  = H2  (V)  + H2  (V)  (1.11) 

t s 

H2  (v) T + — T • (1.12) 

i + (v/v  r i + (v/vc> 

L S 

Example 

Based  on  daytime  tracking  tests  at  WSMR  during  the  CLGP  OT  1 Tests, 
the  following  parameter  values  were  selected  from  Run  75A  for  the  azimuthal 
error : 

V = 0.7  hz  , V = 3 hz 
t 9 s 

a - 0.230  m , a ■ 0.165  m (at  3 km  range) 
t s 

A = 6.8038  10"2  , B = 5.7773  10_3  (m2/hz) . 

Numerical  results  are  shown  in  Table  1 and  in  Figure  1. 

TABLE  1.  ANALYTIC  ESTIMATES  OF  THE  TRACKING  AND  SCINTILLATION 
COMPONENTS  OF  AZIMUTHAL  LASER  SPOT  MOTION  IN  THE  AUTOSPECTRUM 


V 

H2(V) 

/■— s 

p 

N-' 

CM  Cfl 

H2  (V) 

(m^/hz) 

0.3 

6.582 

lo"2 

5.720 

io'3 

7.154 

io-2 

0.5 

5.398 

io-2 

5.621 

io“3 

5.961 

io'2 

0.7 

3.402 

10-2 

5.479 

io-3 

3.950 

io-2 

1.0 

1.317 

io-2 

5.200 

io"3 

1.837 

io-2 

1.5 

3.081 

io"3 

4.622 

io-3 

7.702 

io-3 

2.0 

1.006 

io-3 

4.000 

io"3 

5.005 

io“3 

4.0 

6.375 

io-5 

2.080 

io-3 

2.144 

io-3 

An  alternative  notation  for  the  spectral  density  (autospectrum)  is  r^(v)  , 
for  azimuth,  and  T (v) , for  elevation. 
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GLLD 

model 


est.  from  run  75A 


1 : ■ : ■ T7  ‘ - -I 

Parameters  of  the 
Estimated  Autospectrum 

Length  of  Run 
Lag  Window 
Deg,  of  Freedom 
Nyquist  Frequency 
Bartlett  Band  Width 
90/C  Confidence  Interval/T 


xx 


(V) 


imam 


26,87  sec 
3 sec 

26 

7,5  hz 
0,5  hz 
(0.67,1,69 

, L ..  l i ,iu; 


Frequency,  V (hz) 


Figure  1.  Comparison  of  an  Analytic  Model  of  the  Autospectrum 
of  Azimuthal  Laser  Spot  Motion  With  an  Experimental  Estimate 
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CHAPTER  2 


THE  ANALYTICAL  AUTOCOVARIANCE  FUNCTION  FOR 
LASER  SPOT  MOTION  DEVELOPED 
FROM  THE  SPECTRAL  DENSITY 


The  spectral  density  (autospectrum)  of  laser  spot  motion  in,  say, 
the  x-direction  due  to  the  joint  effects  of  tracking  jitter  and  atmospheric 
scintillation  (beam  steering)  is  given  by 


r (v)  = 


XX 


l+(v/vt) 


l+(v/v  )‘ 

s 


» 0 < v < ”, 


(2.1) 


where 


= 2 Jl 

TT  V, 


2 a 


B =* 


TT  V 


with  variance  contributed  by  tracking  and  variance  contributed  by 


scintillation  o . By  definition, 
s 


2 2 
a + a 
t s 


r (V)  dv. 


XX 


(2.2) 


The  natural  corner  frequency  of  the  second-order  tracking  process  is 


V and  that  of  the  first-order  scintillation  process  is 


V 


An  alternative  expression  for  (2.1)  in  terms  of  the  angular  frequency  to 


is  T1  (to). 

XX 


(U)  ■ rxx 


dv 


dw 


or 
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i 


(2.3) 


rix<“>  - 2rr  rxx  ^ 


r'  (co) 

XX 


A co4  ( 2tt ) — 1 B to2  (2tt)  1 

— + 


4 4 

a)  + to 


2,2 

CO  + CO 
s 


Now,  the  spectral  density  is  by  definition  the  Fourier  cosine  trans 
form  of  the  autocovariance  function  y (t) 


xx 


r'  (03) 


XX 


O 

= " f 

* J n 


Y (t)  cos  cot  dt  . 


xx 


In  general,  the  Fourier  transform  pairs  f(t)  and  g (oo)  are  related  by: 

00 

g(U))  = f(t)  cos  (03t)  dt 


f(t)  - 2tt  g (co)  cos  (cot)  dco  . 


'0 


From  (2.5), 


OO 

= 71-2  [ r*  (co)  cos  (cot)  dco 

****’  I n aa 


^he  following  general  results  will  be  used  to  develop  an  expression 
f°r  ^^(O  for  the  spot  motion  process: 


I (ai  + 032)  1 cos  cot  dco  = (2  TT  a^)""1  e~ait 


0 


4 4 —1 

(a0  + co  ) cos  cot  dco  = 


, - -at/  Jl 

(27T)  a~J  e 2 


TT  a2^ 

sln  <4  + ) • 


(2.4) 


(2.5) 


(2.6) 


(2.7) 


(2.8) 


(2.9) 
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Using  the  expression  for  I*  (to)  in  (2.4), 


» _i  -to  t 

T (to)  cos  tot  dco  = B v (2tt)  e S 
xx  s 


, -to  t//2  to  t 

+ A v (2tt)  e sin  (7-  + ~ -)  . 

t 4 /2 


(2.10) 


With  A and  B as  defined  in  (2,1)  and  with  (2.7), 


v m - S?  n2  "Wtt/v^  . 7T  “t*.  2 "“s’1 

Y„„(t)  - /2  at  e sm  (^  + -*=-)  + ag  e s 


XX 


(2.11) 


The  autocorrelation  function  is  defined  as 


Pxx(t)  - Yxx(t)/Yxx(0)  • 


(2.12) 


In  this  case. 


Y (0)  = a 2 + a 2 

XX  t s 


(2,13) 


Example 

Using  parameter  estimates  from  tracking  run  number  75A  of  the  CLGP 
OT  1 Tracking  Tests, 
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V = 0.7  hz 

V - 3 hz 

s 

o = 0.23  m 

t 

a = 0.165  m 
s 

a 2 / (a  2 + a 2 ) = 0.6586 

t L S 

ag2  / (a2  + ag2  ) - 0.3414 

Pxx(t)  = 0.9314  exp  (-3. lit)  sin  (0 

+ 0.3414  exp  (~18.85t). 
This  result  is  compared  with  t 


= 4.398  r/s 

to  = 18.850  r/s 
s 

(at  3 km  range) 


7854  + 3. lit) 

ie  experimental  estimates  in  Figure 


(2.14) 

2. 


i 
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CHAPTER  3 


DIGITAL  COMPUTER  IMPLEMENTATION  FOR  THE 
FIRST-ORDER  COMPONENT  OF  THE  SPOT  MOTION  ERROR 

Marginal  statistics  for  the  first-order  process  are  assumed  gaussian. 

However,  for  brevity  the  word  "gaussian"  is  often  omitted  in  describing 

the  process.  One  implements  a first-order,  continuous  (or,  analog) 

stochastic  process  with  a first-order  digital  filter  having  a discrete,  mean- 

gaussian  white  noise  input.  The  design  *f  the  filter  is  based  upon  the 

assumption  that  the  continuous-time  process  is  being  time  sampled. 

Continuous-time  (Analog)  Processes 

The  following  first-order  transfer  function  H'(s): 


H'(s)  = 


(B/2ir) 1/2  U)c 
s + 00 


(3.1) 


having  a continuous,  white  noise  input  generates  the  autospectrum 

9 (B/2TT)  U)2 

[H'(C0)]Z 5 ^ , 

or  + w 


(3.2) 


as  in  equation  (2.4)  of  Chapter  2,  where  the  autospectrum  of  the  output  is 
just  the  squared  modulus: 

[H'(to)]2  = H'(jto)  H'(-jw)  (3.3) 

for  a white  noise  input  of  unit  variance. 

For  convenience  a normalized  version  of  (3.1)  is  employed  such  that 
gain  is  unity  for  s equal  to  zero : 


H(s)  = cos/(s  + tOg)  . (3.4) 

Attention  to  the  variance  of  the  input  white  noise  is  deferred  until  later. 


zero 
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Filter  Design 


A time  sampling  of  the  continuous  process  characterized  by  (3.4)  can 
be  used  as  the  basis  of  the  digital  implementation . The  z-transform 
describing  the  digital  implementation  of  (3.4)  can  be  obtained  by  performing 
a bilinear  transformation  from  the  s-plane  into  the  z-plane.  (See  Walsh 
[8]  or  Stanley  [9]).  The  bilinear  z-transform  is  taken  by  substituting 

s = (z  - l)/(z  + 1)  (3.5) 


into  (3.4). 


To  minimize  distortion  of  the  spectrum  of  the  digital  process  relative 

to  the  equivalent  analog  process,  we  use  a warped  analog  cutoff  frequency, 

a (instead  of  u)  ) , where 
s 

a = tan(wgT/2)  = tan(ir  VgT),*  (3.6) 


with  sampling  interval  T and  desired  angular  cutoff  frequency  to  or 
natural  cutoff  frequency  Vg. 


Using  the  warped  cutoff  frequency,  the  digital  transfer  function**  is 

(3.7) 


H (z)  = 

2 Z - 1 

Z + 1 


+ a 


*This  is  a normalized  version  of  the  warping  transformation  of  frequency 
required  by  the  bilinear  transformation.  Normalization  is  achieved  by 
division  by  2/T. 

**The  digital  transfer  function  is  the  ratio  of  the  z-transforms  of  output 
to  input  of  a digital  filter. 

[8]  Walsh,  P.  J.,  A Study  of  Digital  Filters,  AD710381,  Naval  Postgraduate 
School,  Monterey,  CA. , Dec  1969. 

[9]  Stanley,  W.  D.  Digital  Signal  Processing,  Reston  Pub.  Co..  Inc.. 
Reston,  VA. , c.  1975. 
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The  useful  (cannonical)  form  of  (3.7)  Is,  after  manipulation. 


Hz(z) 


(3.8a) 


with 


a 


o 


= a/ (a  + 1) 


b = (a  - 1)/  (a  + 1) . 


(3.8b) 


Remembering  that  z is  the  backspace  operator,  implementation  of  the 


desired  digital  filter  follows  directly  from  (3.8).  Notationally , let 
the  noise  input  to  the  filter  at  time  index  t be  n^_  and  the  output  at 
t be  x , with  t integer. 

Then , 


Ratio  of  Output  to  Input  Variances 

In  selecting  the  value  of  noise  variance  of  the  input  to  yield  a 

particular  output  variance,  it  is  necessary  to  know  the  ratio  of  output 

to  input  variance  of  the  filter.  In  the  following  development  an  expression 

for  this  ratio  will  be  derived.  It  will  be  useful  to  employ  the  following 


(3.9) 


or 


•(3.10) 


Since  the  desired  output  spectrum  requires  a white  noise  input , 


E*-ntnt-k-*  = 0 ’ k ^ 0 


(3.11) 
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notation  for  the  covariance  function: 


Y (k)  = E [x  x 1 ] . 
xx  t t-k 

In  this  notation,  we  desire 


(3.12) 


Y (0 )/a  or  Y (0)/Y  (0) . 

xx  n xx  'nn 


Taking  the  mathematical  expectation  of  the  product  of  n^  with  both 
sides  of  (3.10)  yields 


2 

E[xtnt]  = aQan  > (t  integer),  (3.13) 

since 

E[n  n ] ■ 0 
t t-1 

and 

E[ntxt_1]  = o . 

Similarly,  taking  the  expectation  of  the  product  of  nfc  ^ with  both  sides 
of  (3:10)  yields 


E[xtnt_1]  = aQ(l  - b)  . (3.14) 

After  squaring  both  sides  of  (3.10)  and  taking  the  mathematical  expectation, 


Y (0)  = E[xh  = E[a^  n2  + 2 a2  nn 
xx  t o t o t t-i 

2 2 .,22, 

+ 3o  nt-l  + b xt-l]  ‘ 


2 a bn  x.  - 2 a bn  ,x  , 
o t t-1  o t-1  t-1 


(3.15) 


With  (3.13)  and  (3.14), 


YXX(t»/\n(0>  ‘ 2 Va  + b)- 


(3.16) 
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Alternatively,  with  (3.8b) 


9 


Yxx(0)/Ynn(0)  = a/(a  + 1)  • (3.17) 

This  is  the  desired  result. 

Summary 

In  summary,  to  simulate  a first-order,  mean~zero  process  x^_  having  a desired 
cutoff  frequency  Vg  and  standard  deviation  a , one  can  employ  the  digital 
filter  given  by  (3.10)  with  a gaussian,  white  noise  input  n^_  having  standard 
deviation  given  by 

°n  ' ^nn(0)/YXX(0>l1/2°s 


an  = [(a  + l)/a]1^2as  , 


(3.18) 


with  a given  by  (3.6)  and  filter  constants  given  by  (3.8b). 
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CHAPTER  4 


DIGITAL  COMPUTER  IMPLEMENTATION  FOR  THE 
SECOND-ORDER  COMPONENT  OF  THE  SPOT  MOTION  ERROR 


Analysis  of  a variety  of  man-operated  trackers,  such  as  those  reported 
in  [3],  indicates  that  the  marginal  probability  distribution  function  for 
the  tracking  component  of  laser  spot  motion  is  adequately  described  as 
gaussian.  Additional  experience  with  tracking  records  of  this  type 
indicates  that  a second— order  dynamical  system  characterizes  the  tracking 


error.  For  most  human  trackers  a good  digital  simulation  of  the  mean-zero 
portion  of  the  tracking  error  is  obtained  by  filtering  gaussian  white 
noise  with  a second-order,  low-pass  Butterworth  filter.  As  in  Chapter  3, 
the  design  of  the  digital  filter  is  based  upon  the  assumption  that  the 
continuous-time  (analog)  process  is  being  time  sampled. 

Continuous— time  (Analog)  Processes 

The  following  transfer  function  describes  a second— order  analog  process 

with  analog  corner  frequency  03  and  damping  constant  £: 

a 


0) 


Vs)  - - 


S + 2£  03  S + 03 
a a 


(4.1) 


Specifically  for  low-pass  Butterworth  filters  C = 1 //2  and 


to 

Vs)  = — — 2 — 

s + *^2  to  s + to 


(4.2) 


A 

The  digital  transfer  function  associated  with  (4.2)  is  created  by 
mapping  from  the  s-plane  into  the  z-plane  using  the  bilinear  z-transform. 


The  digital  transfer  function  is  the  ratio  of  the  z-transforms  of  output 
to  input  sequences  of  a digital  filter. 
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(See  Walsh,  P.  J.,  Op  Cit  [8]).  The  bilinaar  z-transform  is  taken  by 
substituting 

s = (z  - l)/(z  + 1) 
in  (4.2).  Thus, 

' Ht  <7TT> 


H (z)  = 


2 2 
0)  (z  + 1) 

a 


(z  - 1)2  + yfl  0)  (z  - l)(z  + 1)  + U)  (z  + 1) 

a • 


(4.3) 


(4.4) 


Reduction  of  H (z)  to  the  form 
z 


-1  j_  -2 

a + a. z + a„z 
o 1 2 


-1  -2 

1 + b^z  + b^z 


(4.5) 


implies* the  digital  filter,  since  z 1 is  equivalent  to  a unit  backspace 
operator.  The  coefficients  in  (4.5)  are  given  by 

•V  ' 2 , 

a = a = U)  /D 

0 z * 

a,  = 2 a 

1 o 

b,  = 2(u2  - 1)/D 

1 a 

b = (U)2  - /T  U)  + 1)/D 

2 a a 


with 


n 

D = 0)  + /2~  0)  + 1 . 

a a 


(4.6) 


Notationally , let  the  noise  input  to  tke  filter  at  time  index  t be  n^ 


and  the  output  at  t be  xt,  with  t an  imte*er.  Then,  from  (4.5), 

*t  + blxt-l  + Vt-2  ■ Vt  + al"t-l  + a2nt-2  • 


(4.7) 


34 


Since  the  analog  angular  frequency  w associated  with  the  bilinear 

ci 

transformation  is  "warped"  or  distorted  relative  to  the  desired  digital 

cutoff  angular  frequency  we  employ  the  relationship  between  these 
parameters: 

w = tan(w  T/2) 
a t 

or 

W = tan  (TTV  T)  , 
a t 

where  T is  the  sampling  period. 

Parenthetically,  note  that  a change  in  the  sampling  period,  for  a 
fixed  digital  cutoff  frequency,  would  change  w and,  by  (4.6),  would 
change  the  coefficients  of  the  digital  filter. 

In  operation,  the  filter  arithmetic  is  performed  as  follows: 


x = a n + a.n  - + a„n  - b x _ h v 
t o t 1 t-1  2 t-2  1 t-1  b2  t-2  * 


with  t = 3,  4,  . . . 


(4.9) 


The  filter  is  initialized  by  assigning  ^ = *2  = 0,  their  mean  values;  and, 
then,  cycling  through  a sufficient  set  of  inputs  to  remove  the  effect  of 
the  initialization  transient. 


The  form  of  (4.9)  is  referred  to  in  the  statistical  literature, 
e.g.,  [10]  and  [11]*,  as  a mixed  autoregressive,  moving  average  model  since 
the  output  xfc  depends  upon  past  values  of  the  output  — xfc  1 and  xfc  2 as 


[10]  Box,  G.  E.  P.  and  Jenkins,  G.  M.,  Time  Series  Analysis:  Forecasting 

and  Control,  Holden-Day,  San  Francisco,  c.  1970. 

Jenkins,  G.  M.  and  Watts,  D.  G. , Spectral  Analysis  And  Its  Applications. 
Holden-Day,  San  Francisco,  c.  1969.  ~ “ 
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well  as  upon  (moving  average)  terms  in  the  input  — nfc  and  nt 

For  the  required  autospectrum  of  nfc  must  be  white.  Thus 


E[ntnt_kJ  » 0,  k 4 0 


a2,  k = 0 
n 


(4.10) 


Furthermore,  future  values  of  the  input  are  uncorrelated  with  the  present 
value  of  the  output,  i*e., 


Ratio  of  Output  to  Input  Variances. 

In  selecting  the  value  of  noise  (input)  variance  to  yield  a particular 
output  variance,  it  is  necessary  to  know  the  ratio  of  output  to  input 
variances  of  the  filter.  An  expression  for  this  ratio  will  be  derived. 

We  employ  a notation  used  by  Jenkins,  Op . Cit. , [11]  for,  respectively, 
the  autocovariance  and  crosscovariance  functions: 


This  result  is  obtained  from  (4.7)  by  successive  multiplication  of  both 
sides  by  nt>  and  and,  then,  by  taking  the  mathematical 

expectation  on  both  sides  in  all  equations.  This  produces: 


E[nt+kxt]  = °»  k>0 


(4.11) 


(4.12) 


In  this  notation,  one  desires 


and 
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Thus, 


Yyn(l)  = (a  - b a )a  . 
xn  1 Ion 


V**(0)  + (1  + b2>  \x(1)  • 

[Vl  + a2(al  - Vo’1  °n  • 

After  squaring  both  sides  of  (4.9)  and  taking  expected  values,  one 
obtains 

(1  - bl  - b2>  \X(0>  * 2 blb2  * 

°n[~2aoalbl  ' 2aoa2b2  " 2a2bl(al  ' blao)  + ao  + al  + a2)- 
One  can  solve  (4.13)  and  (4.14)  for  the  desired  variance  ratio. 

C1B2  ~ BlS 

Y (0)/Y  (0)  = — ■ =—=■ 

Txxv  Tnnvw  A1B2  - BA  » 


(4 


with 


Al  = bl 

a2  = 1 - bj  - b2 


B1  " 1 + b2 


B2  = -2  blb2 


(4 


(4 


C1  = Vl  + a2(al  ~ blao} 

C2  = ~2aoaibl  ~ 2V2b2  " 2a2bl(ai  - blao)  + ao  + ai  + V (4< 
For  numerical  stability  the  denominator  of  the  r.h.s.  of  (4.15a)  must, 
of  course,  be  non-zero.  This  implies  that  the  parameters  b^  and  b2  must 


.13) 


.14) 


.15a) 


15b) 
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lie  within  the  region  defined  by: 


h>2l  <i 


b 


1 


<1 


and 


bx  + b2  > -1. 


This  requirement  is  equivalent  to  the  statement  that  the  poles  of  (4.5) 
lie  outside  the  unit  circle.  See  Box,  Op  Cit , [10],  In  practice,  stability 
requirements  are  satisfied  for  the  problem  of  modeling  human  tracker  error 
using  sampling  rates  of  10  or  20  hertz. 

Summary . 

To  simulate  the  second-order,  mean-zero  component  of  a stochastic 
process  having  the  desired  cutoff  frequency  and  standard  deviation  a 
one  can  employ  the  digital  filter  given  by  (4.9)  with  a mean— zero,  gaussian 
white-noise  input  nt  having  a standard  deviation  given  by 


,1/2 


°n-  ’ <4-«> 

with  variance  ratio  given  by  (4.15)  and  with  filter  coefficients  given 
by  (4.6). 
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CHAPTER  5 


AUTOSPECTRA  OF  THE  DIGITAL  ERROR  PROCESSES 


First  Order  Component 


From  the  digital  transfer  function  of  the  first-order  process  (equation 
(3.8a)), 


with  V defined  over  positive  frequencies  from  zero  to  the  Nyquist  frequency 
Vf*l/(2T). 

The  results  in  (5.1)  and  (5.3)  were  derived  from  a normalized  analog 

transfer  function  having  unity  gain  at  frequency  zero.  Consequently, 

2 

H^(0)  = 1.  Then,  the  autospectrum  of  the  first-order  digital  component, 
r (V),  is  given  by 


where  T (v)  is  the  autospectrum  of  the  digital  white  noise  (input)  process, 
given  by 


H (z)  = 
z 


a (1+z-1) 
o 


(5.1) 


1 + bz 


and  with  the  following  mapping  into  the  frequency  domain: 


jtoT  j 2ttvT 
z = e = e 


(5.2) 


the  squared  modulus 


H?(v)  - H (e^2™1)  H (e-J2”T) 
1 z z 


is 


VV)  = —^2 

1 + b + 2b  cos  2ttvT 


2 

2 a (1  + cos  2?rvT) 


(5.3) 


(5.4) 
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(5.5) 


r„n(v)  =vflc,n 


or,  from  (3.18), 


r (V) 

nn 


2T  a (a+1)  /a 
s 


(5.6a) 


or 


r 

nn 


(v) 


= 2T(l+b) 

2 a2 


2 


a 


s 


0 1 v < vf  . 


(5.6b) 


o 

The  noise  spectrum  T is  comparable  to  the  gain  constant  B in  (1.10). 

As  T approaches  zero  I*  approaches  B. 

Example 

It  is  interesting  to  compare  the  value  of  the  first-order  component 

autospectrum  for  the  analog  process  (equation  (1.6))  with  that  for  the 

equivalent  digital  implementation.  To  this  end,  take  the  parameters 

of  the  example  given  previously  in  Chapter  Is  v =3  hz,  O =0.165  m.  Assuming 

s s 

that  the  sampling  rate  (1/T)  of  the  analog  process  is  20  hz, 

T = 0.05  sec,  and,  by  (3.6), 
a = 0.509525. 


From  (3.8b), 

a = 0.337540 
o 

b = -0.324920  . 

Then,  (5.3)  becomes 


2.  . = 0.227866  (1  + cos  0.314159V) 

V ' 1.1055728  - 0.649839  cos  0.314159V 
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This  result  is  compared  with  that  of  the  corresponding  analog  process 
in  Table  2.  Additionally,  squared  modulii  for  other  sampling  rates  are 
displayed  in  Table  2. 


TABLE  2.  COMPARISON  OF  THE  ANALOG  AND  DIGITAL  SQUARED 
MODULII  OF  THE  FIRST-ORDER  ERROR  PROCESS 


V 

(hz) 

hJ(v) 

analog,  B=1 

2 

H^(v),  digital  at  sampling  rate: 

20  hz 

40  hz 

80  hz 

0.0 

1.0000 

1.0000 

1.0000 

1.0000 

0.3 

0.9901 

0.9915 

0.9904 

0.9902 

0.5 

0.9730 

0.9767 

0.9739 

0.9732 

1.0 

0.9000 

0.9119 

0.9030 

0.9007 

1.5 

0.8000 

0.8183 

0.8045 

0.8011 

2.0 

0.6923 

0.7109 

0.6968 

0.6934 

3.0 

0.5000 

0.5000 

0.5000 

0.5000 

4.0 

0.3600 

0.3297 

0.3531 

0.3583 

6.0 

0.2000 

0.1205 

0.1817 

0.1955 

8.0 

0.1233 

0.0267 

0.0984 

0.1171 

10.0 

0.0826 

0.0000 

0.0545 

0.0755 

r (V)  (m2/hz) 
nn 

5.777  10'3 

8.066  10~3 

7.031  10“3 

6.431  10“3 
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As  noted  from  Table  2,  the  squared  modulii  of  the  analog  and  digital 
processes  are  in  exceilent  agreement  over  the  interval  (0  £ V < Vf/2). 

For  frequencies  within  the  interval  (vf/2  < V < V{)  and,  particularly 
at  the  upper  end  of  this  interval,  the  autospectrum  of  the  digital  process 
departs  significantly  from  that  of  the  corresponding  analog.  Since  there 
is  no  spectral  content  for  the  digital  realizations  above  , i.e.,  all 
of  the  variance  of  the  digital  process  must  occur  at  frequencies  below 
, the  values  of  the  digital  autospectrum  exceed  those  of  the  analog 
having  the  same  variance  for  low  frequencies.  Thus,  the  autospectrum 
of  a digital  process  produced  by  a low  sampling  rate  (and  low  V^)  will 
display  substantial  distortion  relative  to  that  of  the  corresponding 
analog  process.  This  point  is  illustrated  in  Figure  3 in  which  different 
digital  autospectra  of  the  first-order  stochastic  process  are  compared 
to  a corresponding  analog  autospectrum. 
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Autospectrum  (m  /hz) 


a * a fi  -7  h # l 


100 


Frequency,  V (hz) 

Figure  3.  Comparison  of  the  First-Order  Autospectra  for  the  Analog  and  Corresponding 
Digital  Processes 
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Of  course,  if  the  contribution  of  the  first-arder  process  to  the 
total  spot  motion  variance  is  small,  this  distortion  may  be  acceptable. 
However,  if  fidelity  to  the  corresponding  analog  process  is  desired, 
a general  rule  might  be  that  the  sampling  rate  used  in  generating  the 
process  digitally  should  be  20  to  30  times  the  corner  frequency. 

Naturally , the  dynamical  response  of  the  system  accepting  this  first-order 
noise  is  also  a consideration  in  selecting  a sampling  generation  rate. 

Generally  the  dynamic  components  of  a signal  such  as  laser  spot  position 
are  viewed  only  at  discrete  points  in  time,  namely  when  the  pulsed  laser 
flashes.  Thus,  any  spot  motion  due  to  a fundamentally  continuous  (analog) 
process  is  inherently  time  sampled.  Further,  this  sampling  interval 
T*  may  be  different  from  that  used  in  the  digital  implementation,  T.  To 
quantitatively  assess  what  effect  the  digital  generation  interval  T has 
on  the  output  series  obtained  by  sampling  at  a time  interval  T*,  one  must 
exploit  some  sampling  theory.  This  theory  is  developed  below  for  an 
ideal  (delta-function)  sampler  and  applied  to  the  first-order  component 
of  the  stochastic  process. 


The  Transfer  Function  of  a Time-Sampled  Process 


Whenever  a stochastic  process  is  sampled  (or  subsampled),  the  autospectrum 
of  the  sampled  process  may  be  substantially  altered  relative  to  the  original 
process.  This  distortion  will  occur  whenever  the  original  process  has 
an  appreciable  variance  (or  power)  invested  in  frequencies  (v)  above  the 
Nyquist  (or  folding)  frequency  of  the  sampler.  If  the  sampler  is  operating 
at  a sampling  rate  (T*)-1,  the  folding  frequency  will  be  Vf  = (2T*)-1 
and  variance  in  the  original  signal  associated  with  V > Vf  will  be  con- 
founded with  the  variance  associated  with  V for  0 <_  V < Vf . In  this 
section  a theorem  for  ideal  (point)  samplers  will  be  applied  to  several 
processes  and  examples  will  be  offered  to  quantify  the  distortion  accompanying 
sampling.  In  doing  this,  the  author  follows  the  notation  and  results  of 
Kuo  and  Kaiser  [ 12  ] , p.  222  ff . 

Notationally , let  f (nT*)  be  the  value  of  the  original  process  evaluated 
at  points  in  time:  t = nT",  with  n an  integer  and  T*  the  constant  sampling 
interval.  Further,  let  F(s)  be  the  Laplace  transform  of  the  original 
process  and  let  F*(s)  be  the  Laplace  transform  of  the  sampled  process,  i.e.,  of 
the  discrete  series  resulting  from  sampling.  Then,  for  zero  initial 
conditions , 

F*(s)  = E“=_00F(  s + jn2TT/T*)  , (5.8) 

with  s = jco  defined  for  -00  < co  < 00  . 


[ 12  ] Kuo,  F.F.  and  Kaiser,  J.F.  Systems  Analysis  by  Digital  Computer, 
John  Wiley  and  Sons,  New  York,  c.  1966. 
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And, 


F*(ju))  = F(jw)  + Z°°  F (j  (uH-od  ))  + E°*  F ( j (to— a)  ))  , 

11  n n— l n 


with 


w = 2irn(T*)-1 
n 


(5.9a) 


(5.9b) 


and 


(5.9c) 


w = 2irv  . 

Since  the  modulus  of  F ( jto)  will  generally  decline  with  0)  for  values 
of  a)  greater  than  some  value,  say,  w*,  it  will  suffice  to  approximate  F* 
with  a finite  (and  reasonably  small)  number  of  terms  in  the  infinite  sums. 
Thus, 

F*(ju)  s F(jw)  + F(j(omon))  + F(j(03-03n)).  (5.10) 

In  evaluating  the  autospectrum  »f  the  sampled  process  it  will  be 
convenient  to  separately  find  the  real  and  imaginary  parts  of  each  of  the 
terms  in  F*(jw). 

Then, 

Re  {F*( jo))  } = Re  {F ( joo)  > 

£ 

+ Z”  Re  {F(j(aHvd  ))} 
n-j.  n 

+ Zn=l  Re  » (5*1;L) 

and  likewise  for  the  imaginary  part  of  F*(jod),  Im  {F*(jto)}. 

t 

Finally,  the  autospectrum  of  the  sampled  process  is  proportional  to 

I F* ( jto)  | 2 = [Re  {F*  ( joi)  } ] 2 + [Im  {F*  ( ju>)  } ] 2 . (5.12) 

As  a particular  application  of  the  above  results,  take  the  first-order 
analog  process,  for  which 


F(s)  =o)  (s  + w ) 
s s 


-1 


(5.13) 
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In  this  case 


.-1 


F (jto)  = b>  ( jco  + 0)  ) ' + 0)  E , [j  (CJ+W  )+03  ] 


-1 


s n=l 

+ to  E ..  [ j (to— u)  ) + u ] . 

s n=l  J n s 


n 


(5.14) 


Separating  F*  into  real  and  imaginary  parts: 
Re  {F* (jto) } = [ (oj/u)s)2  + l]-1 

oo  2 2—1 

+ E . [ (oj+oj  ) /to  +1] 
n=l  n s 


+ E n [ (to— to  )2/to2  + 1]  1 
n=l  n s 


(5.15) 


and 


- Im  {F"(jto)}  = (w/w  ) [ (to/to  )2  + 1]  1 

s s 


oo  , 2 2 — 1 

+ E ( (otKjO  ) /a)  ) / [((jO+CO  ) / 03  +1] 

n=l v ns'  n s 


+ E°°  . (Cto-to  ) /to  ) C(w-w  ) V + 1]  1 
n=l  ns  s s 


(5.16) 


This  case  has  been  evaluated  for  the  following  specific  parameters: 

V * 1,  3 hz  and  T*  = 0.1,  0.05,  0.025  sec.  The  results  are  dispayed  in 
Figures  4 and  5. 

A second  example  of  the  effect  of  sampling  is  the  case  in  which  a 

digital  implementation  of  a first-order  process  is  subsampled,  i.e.,  in  which 

* -1 

the  sampling  frequency  (T  ) is  a submultiple  of  the  digital  generation 

frequency  T ^ • The  digital  process  with  transfer  function  H (z)  in  equation 

z 

(3.8a)  will  be  the  example  taken  here; 


H (z)  = 


aQ(l  + z ^) 

=1 


(5.17) 


1 + bz 

With  the  mapping  into  the  s-plane: 
sT 

z * e 


(5.18) 
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Frequency,  V (hz> 

Figure  4.  Comparison  of  a First-Order  Anal**  Squared  Modulus  With  the  Squared 
Modulus  of  the  Ideal  Time-Sampled  Series  (Parameters:  V « 1 hz,  T*  = 0.1  sec) 
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|F{jw}|2 


Analog 


I-  - Sampled 
'■7  With  T* 


Sampled 

With  T*  - 0*05  see 


Frequency,  V (hz) 

Figure  5.  Comparison  of  a First-Order  Analog  Squared  Modulus  With  the  Squared  Modulus 
of  the  Ideal  Time-Sampled  Series 

(Parameters:  V = 3 hz,  T*  = 0.05,  0.1  sec) 

s 
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9 


(5.20a) 


F(jw) 


»0d  + exP  (“jwT)  ) 
1 + b exp  (-jojT) 


with  (from  (3.8b)) 


aQ  = a/(a  + 1) 

b - (a  - 1)  / (a  + 1) 

a = tan  (ttv  T) 
s 


In  this  case, 


a (1  + b)  (1  + cos  (joT) 
Re  { F ( jto)  } = -2 _ 

1 + b + 2 b cos  (jl)T 

and 


(5.20b) 


(5.21a) 


a (1  - b)  sin  wT 

-Im  {F(jaj)  } =-2 » (5.21b) 

1 + b + 2 b cos  (i)T 

for  it  _<  wT_<  it  ; and  when  wT  is  outside  of  this  region,  F(jw)  is 
identically  zero. 

Obviously,  the  sampling  interval  T*  must  be  an  integral  multiple  of 
the  generation  interval  T: 

T*  = mT  , m = 1,  2,  3 ... 

In  this  example  F*(jw)  is  evaluated  for  v = 3 hz  and  T = 0.025  and 

s 

T = 0.05  sec  with  m = 2.  Results  are  plotted  in  Figures  6 and  7. 

The  autospectra  for  the  time-sampled  process  in  Figure  7 may  be  compared 
with  the  unsampled  autospectra  shown  in  Figure  3.  The  effect  of  sub- 
sampling order  (m)  on  the  squared  modulus  is  shown  in  Table  3. 


Squared  Modulus 


analog,  w/o  sampling 


^ digital p sampling 
r|te  40  bz 
T “ 0*025  sec 


digital 5 sampling 
rate  20  hz 
T*  * 0*05  sec 


Frequency,  V (hz) 


1*0 


oa 


b 

0.01, I 


0.1 


Figure  6.  Comparison  of  Squared  Modulii  for  Subsampled  Digitally  Generated 
First-Order  Processes  Having  Different  Sampling  Rates 
(Parameters:  Vg  = 3 hz,  T = T*/2) 
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Autospectrum 


impmsi 


Analog 


Digital,  T * 0.05  sec 


DiglEal,  T ■ 0*025  sec 


mmm 


Digital,  T - Q. 0125  sec 


Figure  7*  Comparison  of  the  Autospectra  of  Time-Sampled  Processes:  Analog  and 

Digital  Implementations  (Parameters:  = 3 hz,  ag  = 0.165  m,  T*  * 0.05  sec) 


I 
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TABLE  3.  EFFECT  OF  SUBSAMPLING  ON  THE  SQUARED 
MODULUS  OF  A DIGITAL  IMPLEMENTATION  OF  A 
FIRST-ORDER  DYNAMIC  PROCESS  SAMPLED  AT  20  HZ 

Parameters:  v = 3 hz,  T*  = 0.05  sec 

s 


digital  | F ( jco)  | 2 with  T = T*/m 


subsampling  order  m 


(hz) 

1 

2 

3 

4 

0.1 

0.9987 

0.9987 

1.0311 

1.0537 

0.2 

0.9949 

0.9950 

1.0272 

1.0497 

0.3 

0.9887 

0.9887 

1.0208 

1.0432 

0.4 

0.9800 

0.9801 

1.0120 

1.0342 

0.5 

0.9690 

0.9692 

1.0008 

1.0229 

0.6 

0.9560 

0.9562 

0.9875 

1.0093 

0.7 

0.9409 

0.9412 

0.9721 

0.9937 

0.8 

0.9240 

0.9440 

0.9549 

0.9762 

1.0 

0.8855 

0.8861 

0.9157 

0.9363 

1.2 

0.8420 

0.8430 

0.8714 

0.8913 

1.5 

0.7709 

0.7724 

0.7991 

0.8177 

2.0 

0.6475 

0.6500 

0.6737 

0.6901 

2.5 

0.5306 

0.5342 

0.5551 

0.5694 

3.0 

0.4276 

0.4323 

0.4507 

0.4632 

3.5 

0.3406 

0.3463 

0.3625 

0.3735 

4.0 

0.2687 

0.2753 

0.2898 

0.2995 

4.5 

0.2100 

0.2174 

0.2305 

0.2392 

5.0 

0.1624 

0.1704 

0.1824 

0.1902 
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Aside  from  the  distortions  in  the  autospectrum  (and  process  dynamics) 
due  to  digital  sampling,  there  are  distortions  in  the  spectral  estimates 
created  by  the  statistical  techniques  employed  in  analyzing  time  series 
data.  These  distortions  are  discussed  below. 


Effect  of  Lag  Window  in  the  Estimated  Autospectrum 
In  estimating  an  autospectrum  from  time  series  data,  it  is  essential 
to  provide  a means  of  averaging  or  smoothing  spectral  estimates  to  insure 
stochastic  convergence  as  the  length  of  the  series  (or  record)  grows 
indefinately  large.  For  a single  record,  x(t),  the  smoothing  is  efficiently 
done  by  applying  a weight  function  or  "lag  window",  w(t),  to  the  es  ft  mated 
autocovariance  function,  Y^Ct),  before  taking  the  Fourier  transform  to 
form  the  smoothed  autospectrum,  ^(co).  See  Jenkins  and  Watts,  1968,  Op  Cit.  [U], 
Thus,  for  a continuous  record  of  length  t , 


f (go)  = — 

XX  TT 


Y (t)  w(t)  cos  cot  dt. 


(5.22) 


Without  applying  a weight  function  to  the  theoretical  autocovariance, 
and  for  an  essentially  infinite  record,  one  would  obtain  the  theoretical 
(unsmoothed)  autospectrum: 


OO 


rxx(u)  * I j0  (t>  “ dt. 

Parenthetically , it  is  noted  that  the  spectrum  is  often  expressed  in 
of  the  natural  frequency  V rather  than  the  angular  frequency  to.  Then 


(5.23) 

terms 


9 


r'  (v)  - r (w(v))  1^1 

XX  XX  v 'dv1 

rxx  (V)  = 27T  rxx  (27TV)  • (5.24) 
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Altho  smoothing  is  required  to  reduce  the  variance  of  the  estimated  auto- 

A 

spectrum,  it  does  create  a distortion  of  the  estimate  r^(to)  relative  to 
T (to) . It  is  the  purpose  of  the  developments  of  this  section  to  quantify 

XX 

this  distortion. 

Two  weight  functions  are  frequently  used  in  smoothing:  one  due  to 

Bartlett — 

w (t)  = 1 - t/T  , 0 < t < T 
B m — — m 

=0  , t > T , (5.25) 

m 

where  T is  a lag  parameter,  and  one  due  to  Tukey — 
m 

W (t)  = \ (1  + COs(7Tt/T  )),  0 < t < X 
1 A m — — m 

=0  t > x (5.26) 

m 

In  the  following  analysis  the  simpler,  Bartlett  window  is  used.  However, 
the  effect  shown  on  the  estimated  autospectrum  is  representative  of  several 
lag  windows.  Also  in  this  analysis  we  take  the  estimated  autocovariance 
function  to  be  the  theoretical  autocovariance  of  a first-order  analog  process: 
Y (t)  = Y (0)  CXP  (”Xt) • Specifically,  from  the  second  term  of  (2.11), 

XX  XX 

the  autocovariance  of  the  scintillation  process  is 

Y (t)  = a2  (5.27) 

XX  s 

Then,  for  the  Bartlett  weight  f unction^ the  smoothed  autospectrum  is 


f (to)  =^£s 

XX  7T 


0 T 
2 r m 


f- 

cos(o)t)  e s (1  - t/T  ) dt. 
0 m 


(5.28) 


A 9a 

r (to)  - 


0 00  T 

2 r s m 


(5.29) 


From  (1.10),  the  constant  multiplier  on  the  r.h.s.  of  (5.29)  is  B/2tt. 
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2tt  a 

IT  T <w) 

B xx 


0)  T 
r s m 


to  ™x 

c°s(—  x)  e dx 

s 


(0  T 
s m 

-x  ,0) 

10T  jn  xe  cos(—  x)  dx. 
s m J 0 s 


u 

~t~  |f 


(5.30) 


After  some  manipulation, 


2tt 


B xx 


r (to)  = r (to)  + a (to) 


XX 


-to  T 


+ e s m [a(to)  cos  iot  + 3(to)  sin  tor  1 , 

m m 


(5.31a) 


where 


?xx(»>  = + (co/to  J2)  1 


a (to)  = togTm  1 (to2  - to2)  (to2  + to2)-2 


1 - (to/to  ) 
a (to)  = — — 


“sTm  (1  + <“/“s)2>2 


8 (to)  = - 2 to2x  ^-to (to2  + to2)-2 
sms  ' 


8(to)  = -2  (to/tos)(tosTm)~1  (rxx(to))2  . (5.31b) 

Note  that  the  normalized  smoothed  autospectrum  (2tt/B)  T approaches 

the  unsmoothed  function  as  grows  infinite.  Some  numerical  examples 

of  r were  evaluated  for  several  values  of  the  parameters  V (or  to  ) and 

s s 

Tm*  Results  are  displayed  in  Figures  8 and  9.  The  effect  of  smoothing  is 

to  shift  variance  (or  power)  to  lower  frequencies  and  to  introduce  a damped 

oscillatory  term.  These  distortions  are  not  very  significant  for  values 

°f  V x greater  than  about  1.5. 
s m 
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10  ° 


0.1 


0.01 


Frequency,  v (hz) 

A 

Figure  8.  The  Smoothed  Autospectrum,  F (v)  , Due  to  a Bartlett  Lag  Window  Applied 
to  a First-Order  Autocovariance  With  Associated  Theoretical  (Unsmoothed)  Spectrum 

rxx(v),  (Parameters;  va  * 1 hz,  xm  ■<  l sec) 
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Frequency,  V (hz) 


100 


1*0 


0.01  k 

0.1 


o.i : 


Figure  9.  The  Smoothed  Autospectrum,  ? (v) , Due  to  a Bartlett  Lag  Window  Applied  to  a 

First-Order  Autocovariance  With  Associated  Theoretical  (Unsmoothed)  Spectrum  T (v) 

(Parameters:  V = 3 hz,  T =0.5  sec) 

s in 
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Second-Order  Component 

The  digital  transfer  function  used  in  describing  the  second— order , 
tracking  component  of  the  laser  spot  motion  (equation  (4.5))  is  repeated 
here. 

, -1  . -2 

a + a-  z + a9  z 

H (z)  = -2 ~r l—Jr  ‘ (5.32) 

z 1 , 1 1 i i 2 

1 + b^  z + z 

The  autospectrura  of  the  process  is  proportional  to  the  squared  modulus 

|H  (exp (jwT) ) | 2 . Notationally, 
z 

H? (to)  = |H  (eja)T)|2  = |H  (ej27rVT)|  2,  (5.33) 

z z z 

with  sampling  generation  interval  T. 

Using  De  Moivre's  theorem, 

= cos  toT  + j sin  toT,  (5.34) 

with  (5.32)  and  (5.33), 

2 2 
H„  (to)  = [ (a  + a,  cos  toT  + a.  cos  2toT) 

L O 1 L 

+ (a^  sin  toT  + a^  sin  2toT)^']/ 

[(1  + cos  toT  + b?  cos  2toT)2  + (b?  sin  toT  + b2  sin  2toT)2].  (5.35) 


Using  equation  (4.6)  and  after  some  manipulation,  (5.35)  becomes 

A 

9 co  (3  + 4 cos  coT  + cos  2coT) 

»2(“)  ‘ —l, 4 4 ’ <5-36> 

3(oi  + 1)  + 4(0)  - 1)  cos  0)T  + (oi  + 1)  cos  2o)T 

a a a 

where 

co  - tan  (ttv  T).  (5.37) 

a t 

Either  (5.35)  or  (5.36)  can  be  used  to  evaluate  the  squared  modulus  of 
the  second-order  dynamic  component.  Altho  (5.36)  has  a simpler  form 
than  (5.35),  computational  experience  using  single-precision  arithmetic 
on  the  IBM  360  series  computers  has  demonstrated  that  (5.36)  is  somewhat 
more  sensitive  to  truncation  error  than  (5.35).  In  fact,  it  is  noted 
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generally  all  second  order  transfer  functions  and  squared  inodulii 
are  more  sensitive  to  loss  of  precision  due  to  arithmetic  truncation 
than  are  their  first-order  counterparts.  When  using  IBM  360  computers 
It  is  recommended  tnat  double-precision  arithmetic  be  used. 

The  above  expressions  for  the  second— order  squared  modulus  are 
normalized  so  that  H^(0)  = 1.  To  calculate  the  autospectrum,  T (v)  , 
of  the  second-order  (tracking)  process,  one  employs: 

r^Cv)  - h2<2«)  rnn(v)  , (5.38) 


where  ^(v)  is  the  input  noise  spectrum  to  the  second-order  digital 
filter.  The  noise  spectrum  is  a function  of  the  sampling  generation  in- 
terval, T,  and  the  variance  of  the  tracking  process,  a2: 


'nn  - 21  * 

with  the  ratio  of  variances  Y (0)/y  (0) 

nn  xx 


(5.39) 


given  by  (4.15a). 


As  with  the  first-order  process,  it  is  useful  to  make  numerical 
comparisons  between  the  squared  modulus  of  the  digital  implementation 
of  the  second-order  (tracking)  component  and  that  of  the  corresponding 
analog  process.  The  transfer  function  for  a general  second-order  analog 
process,  H^(s),  is  given  by  equation  (4.1),  The  general  squared  modulus 
is,  then, 


H^(jw)  = 1/[(1  - (w/cja)2)2  + 4C2  (w/wa)2] 
With  the  Buttefworth  assumption  of  £ = 1/^2  and  with 

H2(V)  = 1/[1  + (v/vt)4]  , 


to  = 


2TTVt, 


(5.40) 


(5.41) 


as  in  equation  (1.1)  with  A = 1. 

The  analog  gain  constant  A (equation  (1.5))  represents  the  analog 
input  noise  spectrum  and  is  comparable  with  r given  by  (5.39).  In  fact. 


lim  T 
T-K) 


nn 


= A . 


(5.42) 
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Example 


Pursuing  the  example  of  Chapter  1,  with 

V = 0o7  hz 
t 

a = 0.23  m 
t 

T = Ool  sec, 

= 0.219912,  and  from  (5.36): 

H^Cv)  = 2.33879  10'3[3  + 4 cos(0, 62832V) 

+ cos(l. 25664V) ]/ [3.00716 

-3.99064  cos(0. 62832V)  + 1.002339  cos(1.25664v) ] . (5.43) 

Numerical  results  from  this  expression  and  similar  digital  squared  moduli! 
calculated  with  T = 0.05  sec  and  0.025  sec  are  compared  with  the  analog 
squared  modulus  in  Table  4.  Spectral  amplitudes  for  the  noise  are  also 
shown  in  Table  4.  Spectra  are  plotted  in  Figure  10. 


TABLE  4.  COMPARISON  OF  THE  ANALOG  AND  DIGITAL  SQUARED 
MODULI I OF  THE  SECOND-ORDER  (TRACKING)  ERROR  PROCESS 


V 

(hz) 

analog 

H^(V) 

4=1 

digital,  H2(v)  with 

T (sec) : 

0.025 

0.05 

0.10 

0.1 

1.000 

1.000 

1.000 

1.000 

0.2 

0.993 

0.993 

0.993 

0.994 

0.3 

0,967 

0.967 

0,968 

0.969 

0.4 

0.904 

0.904 

0.905 

0.907 

0.5 

0.793 

0.794 

0.794 

0.799 

0.6 

0.649 

0.650 

0.650 

0.653 

0.7 

0.500 

0.500 

0,500 

0.500 

0.8 

0.370 

0.369 

0.368 

0.365 

1.0 

0.194 

0.193 

0.191 

0.183 

1.2 

0.104 

0.103 

0.101 

0.092 

1.5 

0.0453 

0.0447 

0.0428 

0.0357 

2.0 

0.0148 

0.0144 

0.0132 

0.0089 

3.0 

0.00296 

0.00275 

0.00220 

0.00070 

5.0 

0.00038 

0.00031 

0,00015 

0.00000 

rnn0u)  (m2/hz) 

0.0680 

0.0682 

0.0685 

0.0695 
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digital 


digital 


analog 


Autospectrum  (m  /hz) 


itnf! 


«f 


Frequency,  V (hz) 

Figure  10.  Autospectra  of  the  Second-Order  (Tracking)  Process (Parameters : 
V = 0.7  hz,  o = 0.230  m,  T = 0.1  sec) 


CHAPTER  6 


SPECTRAL  MOMENTS  AND  PARAMETER  ESTIMATION 


The  spot  motion  process  has  been  shown  to  consist  of  a mixture  of 

first- and  second— order  dynamic  components  characterized  by  four  parameters — 

2 

the  two  crossover  frequencies  v and  V.  and  the  two  variances  a and 
2 st  s 

°t  or,  alternatively,  the  spectral  gain  constants  B and  A.  The  existence 
of  four  parameters  to  estimate  from  the  autospectrum  suggests  that  ex- 
pressions for  the  first  four  spectral  moments  (zeroth  thru  third)  in 
terms  of  the  four  unknowns  could  be  developed  and  used  with  experimental 
values  of  the  spectral  moments  to  solve  for  the  estimates.  This  concept 
was  explored  and  found  to  be  impractically  complex  for  a mixture  of  dynamic 
components.  However,  for  evaluating  processes  having  only  one  dynamic 
component,  the  concept  of  matching  spectral  moments  to  produce  parameter 
estimates  appears  to  have  promise.  The  estimating  relations  for  a pure 
Butterworth  second— order  process  and  for  a pure  first-order  process 
are  derived  below. 


Spectral  Moments  for  the  Second-Order  Butterworth  Process 

The  autospectrum  for  the  continuous-time  (analog)  second-order 
Butterworth  stochastic  process  {x(t)}  is  given  by 


with 


r (v)  = 

xx 


4 O 


1 + (v/vt) 


, 0 _<  V < °°, 


IT 


V 


(6 . la) 


(6.1b) 


By  definition  the  k th  spectral  moment  is  given  as 


X 


k 


f vk  T (v)  dv 
j o xx 


9 


or 


, * k+l 

X,  = A V 
k t 


k 

x dx 
0 1 + x^ 


(6.2) 
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Using  the  general  result: 


•O 


X 

X 

X 


0 

1 

2 


x*3  ^ dx 
1 + xq 


= ('U'/ q)  tl/sin(p7r/q)  ] , 


1 1 p < q , 


= A vt  (2)  JT 

■ * ”t2  <?> 

= A vt3  (2)  («/?) 


(6.3) 


(6.4) 


For  all  k>3,  the  moments  of  the  analog  process  are  infinite*  However, 
the  parameters  A and  can  be  estimated  using  and  X as  follows: 


A 


V 


t 


2 V 

7 T X^J^ 

yfT  X 


1 

0 


(6.5) 

(6.6) 


In  practice  A and  vt  would  be  developed  by  substituting  the  sample 

A A 9 A 

estimates  XQ  (or  op  and  for  XQ  and  X , respectively. 


The  moments  for  a digital  implementation  of  this  process  are,  of 
course,  finite.  Let  be  the  Nyquist  (folding)  frequency.  Then,  the 
digital  moments  are  approximated  by 


or 


X 


T 

xx 


(v)  dv 


I 


(6.7) 


1 = A V 
k 


k 

k+1  [ xK  dx 

* Jo  1 + x4  ’ 


(6.8) 
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where 


5 = vf/vt  . 


(6.9) 


For  notational  convenience,  define  the  k th  normalized  moment  A,  as 

k 


Ak  = Xk/(A  v*+1)  , k>0. 


(6.10) 


Then, 


Xn  = [*-(*,  + ^ ^ h * 2 tan'1 


0 4 /2 


5 - ^2  5 + 1 


r)] 


1-5 


Y _ 1 „ -l/r2. 

*12  tan  & ) 


- -h  tfa(V-^  £ + *>  ♦ 2 tan-1^)) 


2 4/2 


5 + 5 + l 


l-C 


X3  = \ £n(l  + £4) 


(6.11) 


Spectral  Moments  for  the  First-Order  Process 

The  autospectrum  for  a first-order  process  {x(t)}  with  crossover 

frequency  V and  gain  constant  B is  given  by 
s 


r 

XX 


(V) 


r-  , 0 _<  V < 00  . 

i + (v/v  r 
s 


(6.12) 


From  the  definition  of  spectral  moment  (equation  (6.2))  and  the  general 
result  of  (6.3),  it  is  seen  that  only  the  zeroth  moment  of  the  first- 
order  process  is  finite: 


A 


0 


B dv 

01+  (v/v  )2 
s 


AQ  = B(tt/2) 


(6.13) 
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As  with  the  second-order  process,  define  the  moments  of  the  digital 
implementation  of  the  first-order  process  (Xk>  as  definite  integrals 
with  the  Nyquist  frequency  as  the  upper  limit.  Thus, 


% 

X, 


-C 


v r (v)  dv 

XX 


K = B (V  )k+1 


f xk  dx 
JO  1 + x2  * 


(6.14) 


with 


n - vf/vs 


(6.15) 


The  normalized  moments  for  the  first-order  process  are  defined 


as 


\ = V(B  vs+1)  * k > °* 

For  the  first-order  process, 

XQ  = tan_1(n) 


(6.16) 


\ \ -£n(l  + n2) 

X = D - tan  1(n) 


X3  = 2~  “ 2 + I2) 


(6.17) 


By  inspection  of  (6.16)  and  (6.17)  one  can  write  the  following 

two  equations  which  can  be  solved  simultaneously  for  B and  v : 

s 
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(6.18) 


- 2 <&*,> 

2 [tan  1(n)]Z  0 1 

^ Oj  1 /o 

= vf[(B  vf-XQ)/X 2]lu, 


with 


vs  - vf/n  . 


(6.19) 
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APPENDIX  A 

SUBROUTINE  TO  GENERATE  LASER 
SPOT  MOTION 


Subroutines  SPOTIN  and  SPOTMO— an  entry  point  in  SPOT IN-- generate 
spot  motion  in  azimuth  and  elevation  coordinates  for  use  in  laser  signature 
models.  SPOTIN  initializes  the  program  SPOTMO  by  providing  azimuth 
and  elevation  standard  deviations  for  both  the  tracking  and  scintillation 
stochastic  components  as  well  as  the  crossover  (corner)  frequencies 
of  each  of  the  components.  The  sampling  generation  interval  is  also 
specified.  For  subsampling,  subroutine  SPOTMO  is  called  several  times 
for  each  repetition  of  the  laser  interpulse  interval.  Thus,  the  sampling 
generation  interval  must  be  chosen  in  a manner  appropriate  to  the  use 
of  this  routine.  SPOTIN  computes  the  filter  coefficients  for  the  tracking 
component  by  means  of  equation  (4.6)  and  the  coefficients  for  the  scintillation 
(first-order)  component  via  (3.8b). 

The  input  requirements  of  both  SPOTIN  and  SPOTMO  are  specified 
in  the  following  source  program  listing.  Output  from  SPOTMO  is  also 
specified  in  that  listing.  In  operation  SPOTIN  is  called  only  once  for 
initialization  and  SPOTMO  is  called  where  the  output  position  array 
is  required,  typically  nested  in  a loop  of  the  calling  program. 
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on  non  on  non  ooooooooonoonoo  on 


<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  SPOTIN  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
--  SPOT  MOTION  DRIVER  WITH  SCINTILLATION  >>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
SUBROUTINE  SPOTIN  (SIGMA.  CF . DT) 

INTEGER  TR/I/.  SC/2/ , AZ/1/.  EL/2/ 

dimension  Sigma (2*2)  * cf(2) 

— SPOTIN  INITIALIZES  THE  ROUTINE  SPOTMO  FOR  THE  2-DIMENSIONAL  (AZI- 
MUTH 5.  ELEVATION)  MOTION  OF  THE  DESIGNATOR  BEAM  ABOUT  THE  DESIRED 
DIRECTION. 

INPUT  REQUIREMENTS: 

* SIGMA  --  ARRAY  OF  STANDARD  DEVIATIONS  OF  ERROR  COMPONENTS. 

UNITS  SAME  AS  "OUTPUT"  OF  SPOTMO.  BEL°W. 

SIGMA(I.I)  = SIGMA(AZ.TR)  = S.  D.  IN  AZIM  OF  TRACKING  PROCESS. 

SIGMA (2*1)  = SIGMA  (E.L«TR>  = S.  D.  IN  ELEV  OF  TRACKING  PROCESS. 

SIGMA (1*2)  = SIGMA(AZ.SC)  = S.  D.  IN  AZIM  OF  SCINT  PROCESS. 

S I GMA (2.2)  = SIGMA(EL.SC)  = S.  D.  IN  ELEV  OF  SCiNT  PROCESS. 

* CF  --  VECTOR  OF  CORNER  FREQUENCIES  OF  PROCESSES  (HZ). 

CF ( 1 ) = CF(TR)  = C.  F.  FOR  TRACKING 

CF ( 2 ) = CF(SC)  = C.  F.  FOR  SCINTILLATION 

* DT  --  SAMPLING  INTERVAL  ( SEC  > 

DIMENSION  XT (2*3) , ZT(2,3),  XS(2,2).  ZS<2,2).  OUTPUT(2,3) 

DATA  PI  / 3.141593  /.  NET  /3/ 

--  CHECK  FOR  PRESENCE  OF  TRACKING  ERROR 

ITRACK  = 0 

IF  ( ( (SIGMA(AZ.TR) .LE.O.)  .AND.  (SIGMA(EL.TR)  .LE.O.)  ) .OR. 

1 (CF(IR) .LE.O.)  ) GO  TO  150 

ITRACK  = I 

--  TRACKING  ERROR  IS  present 

— CLEAR  TRACKING  NOISE  ARRAYS  AND  COMPUTE  FILTER  COEFFICIENTS 

DO  100  1=  AZ.  EL 
DO  100  J=  1 . 3 
XT ( I , J)  = 0. 

100  ZT ( I . J)  = 0. 

ZETAAL  = SORT (0.5) 

OMEGAA  = TAN (PI«CF (TR) *DT) 

DUMMY  = I.  / ( 1 . ♦ OmEGAA* (2.  *ZETAAL  ♦ OMEGAA)  ) 

A0  = 0MEGAA**2  * DUMMY 
Al  = A0  ♦ A0 

BI  = 2.  * (OMEGAA**2  - 1.  ) * DUMMY 
B2  = (I.  - OMEGAA* ( 2 • *ZETAAL  - OMEGAA)  ) * DUMMY 
SDNT  = SORT  ( (B2-l.)*( (H2+I.)**2  - B1*9I)  / < A0*A 1 * < ( B2- 1 . ) **2  - 
1 ( B 1 -2 . ) **2 ) ) ) 

TMAX  = 5 . *SDNT 
TM IN  = -TMAX 
150  CONTINUE 

--  CHECK  FOR  PRESENCE  OF  SCINTILLATION 
ISCINT  = 0 

IF  ( ( (SIGMA (AZ. SC) .LE.O.) .AND. ( S IGMA ( EL . SC) .LE.O.) ) .OR. 

I (CF (SC) .LE.O. ) ) GO  TO  250 

ISCINT  = I 

--  SCINTILLATION  IS  PRESENT 

--  CLEAR  SCINTILLATION  NOISE  ARRAYS  AND  COMPUTE  FILTER  COEFFICIENTS 
DO  200  I = AZ,  EL 
DO  200  J = 1.2 
XS( I , J)  = 0. 

200  ZS(I.J)  = 0. 


70 


non  o non  oooooooooooooo 


OMEGA A = TAN(P1*CF(SC)*DT) 


AAO  = OMEGAA  / ( OMEGAA* 1 * ) 

BB  = (OMEGAA-l.)  / (OMEGAA+1.) 
SDNS  = SORT ( 1 « + 1. /OMEGAA) 
SMAX  = 5,*SDNS 
SMIN  = -SMAX 
£50  CONTINUE 
RETURN 


ENTRY  SPOTMO  (NFWD*  1SEEDT,  1SEEDS,  OUTPUT) 

1001  FORMAT  ( 'O^oERROR  --  ',  1 1 0 * *-TlME-STEP  ADVANCE  CALLED  FOR  IN  SPO 
1 TMO  — RUN  ABORTED. • ) 

— SPOTMO  COMPUTES  SPOT  MOTION  COMPONENTS. 

INPUT  REQUIREMENTS: 

* NFWD  — NUMBER  OF  TIME-STEPS  FORWARD  TO  ADVANCE. 

* ISEEDT  — RANDOM  NUMBER  SEED  FOR  TRACKING 

* 1 SEEDS  — RANDOM  NUMBER  SEED  FOR  SCINTILLATION 

* OUTPUT  — ARRAY  OF  OUTPUT  SPOT  MOTION  COMPONENTS 

OUTPUT  (1*1)  = OUTPUT  (AZ*TR)  = AZ1M  ERROR  FROM  J.RACK1NG  PROCESS 
OUTPUT (2*1)  = OUTPUT (EL*TR) 

OUTPUT (1*2)  = OUTPUT (AZ, SC) 

OUT  PUT ( 2 * £ ) = OUTPUT (EL*SC) 

OUTPUT (1*3)  = OUTPUT (AZ*NET)  = SUM  OF  TRACK  6 SCINTILLATION 
OUTPUT (2*3)  = OUTPUT (EL, NET) 


- CHECK  ADVANCE 

IF  (NFWD.GT.O)  GO  TO  300 
WRITE  (6*1001)  NFWD 
STOP 

300  CONTINUE 


CHECK  FOR  TRACKING 


IF  (1TRACK.EQ.1)  GO  TO  350 
OUTPUT (AZ*TR)  = 0. 

OUTPUT (EL*TR)  = 0. 

GO  TO  400 
350  CONTINUE 

DO  380  K=1,NFWD 

- BACKSPACE  NOISE  ARRAYS  AND  GENERATE  NEW  NOISE 

DO  370  1 = AZ*  EL 
DO  360  J=l,2 
JJ  = 4-J 

XT ( 1 * JJ)  = XT ( 1 * JJ-1 ) 

360  ZT ( I * JJ)  = ZT ( 1 , J J- 1 ) 

CALL  NORMXX  ( TM1N  * TMAX  , 0.*  SDNT  * ZT(1*1),  13EEDT) 

370  XT  (1,1)  = A0*(ZT < 1,1 ) ♦ZT ( 1,3) ) ♦ Al*ZT < 1 *2) -B1*XT < I »2> - B2*XT<1,3) 
380  CONTINUE 

OUTPUT (AZ*TR)  = XT ( AZ , 1 ) *S I GMa ( AZ * TR ) 

OUTPUT (EL, TR)  = XT (EL* 1 ) #S1GMA (EL* TR) 

400  CONTINUE 

- CHECK  FOR  SCINTILLATION 

IF  ( ISCINToEQ. 1 ) GO  TO  450 
OUTPUT (AZ, SC)  = 0. 

OUTPUT (EL, SC)  = 0. 

GO  TO  500 
450  CONTINUE 

DO  480  K= 1 *NFWQ 
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non  no  non  on  non 


c— 


470 

480 


500 


BACKSPACE  NOISE  ARRAYS  AND  GENERATE  NEW  NOISE 
DO  470  1=  AZ,  EL 
XS ( I * 2)  = XS ( I * 1 ) 

ZS ( I * 2 ) = ZS ( I * 1 ) 

CALL  NORMXX  (SMIN  , SMAX  , 

XS(I.I)  = AA0*(ZS(T.I>*ZS(I»2) ) 

CONTINUE 

OUTPUT ( AZ  * SC ) = XS(AZ,I)"SIGMA(AZ,SC) 

OUTPUT (EL. SC) 

CONTINUE 


o..  sdns,  zs(i,d  , i seeds > 

- BB^XS (1,2) 


= XS (EL , 1 ) #S 1GMA (EL « SC ) 


DEFINE  NET  SPOT  MOTION 


OUTPUT (AZ, NET)  = OUTPUT ( AZ , TR ) 
OUTPUT (EL, NET)  = OUTPUT ( EL , TR ) 
RETURN 
END 


OUTPUT (AZ, SC) 
OUTPUT (EL, SC) 


>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  NORMXX  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 

NORMXX  --  NORMXX  --  NORMXX  --  NORMXX  — NORMXX  — NORMXX  00249000 


SUBROUTINE  NORMXX  (AMIN,  AMAX,  AMEAN,  SIGMA,  X,  I SEED)  00249100 

00249200 

THIS  SUBROUTINE  GENERATES  A NORMAL  DEVIATE  00249300 

00249400 

CALL  RANDMM ( I SEED , X ) 00249500 

2 = X 00249600 

CALL  RANDMM ( ISEED.X)  00249700 

X = ( ( (-2.0"ALOG(Z>  ) "<>0.5)  <MCOS  (6.283*X) ) )*SIGMA  ♦ AMEAN  00249800 

IF  (X.LT.AMIN)  X = AMIN  00249900 

IF  (X.GT.AMAX)  X = AMAX  00250000 

RETURN  00250100 

E-ND  00250200 


>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  RANDMM 
--  RANDMM  — RANDMM  — RANDMM  — 
SUBROUTINE  RANDMM ( I SEED , X ) 

this  subroutine  generates  uniform 

I SEED  = ISEED*65539 
IF ( I SEED) 2266,2266,2277 
2266  ISfED  = I SEED  ♦ 2147483647  ♦ I 
2277  X = I SEED 

X = X* . 46566 1 3E-9 

RETURN 

END 


<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 

RANDMM  — RANDMM  — RANDMM  00250300 

00250400 

00250500 

DEVIATES  00250600 

00250700 

00250800 

00250900 

00251000 

00251100 

00251200 

00251300 

00251400 


72 


APPENDIX  B 
PROGRAM  TO  EVALUATE 
THE  AUTOSPECTRUM  OF  TIME-SAMPLED 
FIRST-ORDER  PROCESSES 


Purpose 

This  is  a utility  program  which  evaluates  squared  modulii  associated 
with  time-sampled  first-order  analog  and  digital  stochastic  processes. 

The  phase  of  the  time-sampled  digital  process  is  also  calculated.  Out- 
puts are  calculated  for  a prescribed  set  of  input  frequencies. 

Description 

The  squared  modulus  of  a time-sampled,  digitally-generated,  first-order 
stochastic  process  is  calculated  by  means  of  equations  (5.10),  (5.12), 

(5.21a),  and  (5.21b).  The  infinite  series  is  truncated  at  n*,  as  in 
(5.10).  The  value  of  n*  is  determined  implicitly  by  requiring  that 
the  largest  relative  error  of  the  real  and  imaginary  components  of  the 
n term  in  the  series  be  smaller  than  a specified  (input)  relative  error. 

For  comparison,  the  autospectrum  of  an  unsampled  analog  process 
having  the  same  specified  corner  frequency  is  also  calculated  and  printed. 
Additionally,  the  squared  modulus  of  the  tion-subsampled  digital  implementation 
is  calculated  and  printed.  The  latter  derives  from  the  F(jco)  term  in 
equation  (5.10).  The  phase  angle,  <j> , of  the  time-sampled  digital  process 
is  calculated  from 

<t>  " tan  1[Im{F*(jco)}/Re{F*(jw)}] . 


Input 


The  program  requires:  (a)  a descriptive  alphameric  title  card; 

(b)  the  desired  process  standard  deviation  (used  in  calculating  the  analog 
autospectrum);  (c)  the  sampling  time  interval;  (d)  corner  frequency; 

(e)  a maximum  relative  error  in  the  subsampled  digital  squared  modulus; 

(f)  the  number  of  frequencies  for  which  results  are  to  be  calculated; 

(g)  the  order  of  subsampling  or  ratio  of  sampling  interval  to  digital 
generation  interval;  and  (h)  the  array  of  discrete  frequencies  for  which 
results  are  to  be  calculated. 

Output 

Examples  of  program  output  are  shown  following  the  source  program 
listing.  Output  titles  are  generally  self-descriptive.  The  input  data 
are  immediately  echoed.  For  clarification  note  that  SQ  MODUL  0 refers 
to  the  squared  modulus  of  the  non-subsampled  digital  process.  Except 
for  truncation  error,  the  results  of  SQ  MODUL  0 are  identical  to  SQ  MODULUS 
whenever  the  subsampling  order  is  unity. 
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31 
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34 

35 


3>  JOB 

C — 

C 

c * # * «■ 
c **#  «■ 
q # « * # 

q -turn# 

Q«il 
Q # <MH» 

c * * * « 

C * * * * 
C«««  « 

Q #**■» 
C «■#«■* 

c «■*«■•» 


11 

100 

200 

300 

400 

1 

500 
c •» 

600 

1 

(,  * # * « 


Q tt  tt  Ott 


Q tt 


RnH,KP=29»RUN=FREE»LINES=5* 


EVAL  OF  SAMPLE*  TRANSFER  FUNCTION  

PROGRAM  TO  EVALUATE  THE  CHARACTERISTICS  OF  A 
TIME-SAMPLED  FIRST-OPDER  DIGITAL  DYNAMIC  PROCESS. 

Inputs: 

TITLE  = AN  ALPHAMERIC  DESCRIPTION  OF.  THE  FUNCTION 

STIJV  = STANDARD  DEVIATION  OF  THE  DIGITAL  PROCESS  BEFORE  SAMPLING 

T = SAMPLING  PERI  OP  (INTERVAL) 

FS  = CORNER  FREQUENCY 

EPS  = MAX  RELATIVE  ERROR  IN  THE  SQUARED  MODULUS 
MMX  = THE  NUMBER  OF  ENTRIES  IN  THE  FREQUENCY  ARRAY 
MULTIP  = THE  ORDER  NUMBER  OF  SUBSAMPLING 
= 1*  2*  3 • ETC. 

FREQ(M)  = THE  ARRAY  OF  DISCRETE  FREQUENCIES  AT  WICH  THE 
TRANSFER  FUNCTION  IS  TO  BE  EVALUATED. 

THE  ORIGINAL  (ANALOG)  SQUARED  MODULUS  IS  ALSO  CALCULATED. 

I GEN  = SAMPLING  GENERATION  INTERVAL  FOR  THE  DIGITAL  PpOCESS*  T GEN: 
DIMENSION  FREQ<50) ,TITLE(20) 

DATA  TWOPI/6, 283185/ 

CONTINUE 


t/multip 


READ  (5*100 * END  = 1 ) TITLE 
FORMAT (20A4) 

WRITE  (6,200)  TITLE 
FORMAT ( 1H1 ,20A4) 

READ  (5,300)  STDV,T,FS, EPS, MMX, MULTIP 
FORMAT (4F 10-0,213) 

WRITE  (6,400)  STDV,T,FS, EPS, MMX, MULTIP 

FORMAT ( 1H0,6X,4hSTDV,9X, 1HT ,8X,2HFS,7X,3hEPS,2X,8HMX  FREQS, 
10H  ORU  S SMP/1H  ,3F10.4,F10.6,7X, I3,7X, 13) 

READ  (5,500)  (FREQ(I) ,1=1 ,MMX) 

FORMAT (8F10.0) 

WRITE  HEADINGS 
WRITE  (6,600) 

FORMAT  ( 1H0,6X,9HF'REQ  ( HZ ) , 3X  , 1 2H  ANAL  AUTO^PK  ,5X  , 1 OHSQ  MODUL 
5X  » 1 OHSQ  MODULUS, 8X , 7HMODULUS » 6X , 9HPHASE  ( D ) ) 
OMEGS=TWOPl#FS 

COMPUTE  THE  CONSTANTS  FOR  THE  DIGITAL  TRANSFER  FUNCTION 
TGEN=T/FLOAT (MULTIP) 

PI=TWOPI/2. 

AF  =ATAN  <PI*FS*TGEN) 

AO  = AF / ( AF+  1 . ) 

8=(AF-1,)/(AF*1.) 

CONREL=AO*(l.+B) 

CONIMG=-AO# ( 1 .-B) 

CONC=2.»B 
COND=l  ,+B<n>2 

GA I N=SQRT (?.»TGEN/AO) »STDV 
start  frequency  LOOP 
DO  10  M* 1, MMX 
F = F REQ ( M ) 

OMEG=TWOPI*F 
COSX=COS (OMEG*TGEN) 

SINX=SIN (OmEG*TGEN) 

IPASS=0 

COMPUTE  THF  REAL  PART  OF  THE  T RANSFER- FUNCT I ON 

DENOM=COND+CONC*COSX 

SUMREL=CONREL<* ( 1 • +COSX ) /DENOM 

savei=sumrf;l  • 


0* 
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SUMIMG=CONIMG*SINX/UENOM 

SAVE2=SUMIMG 

accumulate  OVER  POSITIVE  and  negative  harmonics 

DO  20  N=1 .24 

IS-WPOS  = 0 
I SWNEG=0 

OMt  GN  = TWOP  I OFLOAT  <M)/T 
TERMP=0. 

TEKMN=0 . 

OMEGPT= (OMFG+OMEGN) oTGEN 
OMEGMT  = (OMEG-OMEGN) #TGEN 
IF (OHEGPT.GT.PI ) I5WP0S=1 

IF  (OMEGMT  .I.T.-PI  .OR. OMEGMT. GT. PI  ) ISWNEG=1 
IFUSWPOS.EO.l)  GO  TO  30 
COSP=COS {OMEGPT) 

DENOMP=COND+CONC«COSP 
TERMP=CONRFL# ( 1 . + COSP ) /DENOMP 
30  CONTINUE 

IF ( ISWNEG.EO. 1 ) GO  TO  35 
COSN=COS (OvEGMT ) 

DEMOMN=C()Nf)*CONCttCOSM 
TEPMN-CONRFL* < 1 • ♦COSN) /DENOMN 
35  CONTINUE 

R=TERMP+TEPMM 

SUMREL=SUMREL+P 

IF (R/SUMREL.LT.EPS.AND.N.GT.4)  IPASS=1 

COMPUTE  THF  NEGATIVE  OF  THE  IMAGINARY'  PART  OF  THE  TRANSFER  FUNCTION 
TERMP=0. 

TERMN=0. 

IF (1SWPOS.EO. 1 ) GO  TO  40 
TERMP  = CONIMG*SIN  (OMEGPT) /DENOMP 
40  CONTINUE 

IF ( ISWNEG.EO. 1 > GO  TO  45 
T ERMN=C0NIKG*S1N (OMEGMT ) /DENOMN 
45  CONTINUE 

R=TERmP+TERMN 

SUM1MG=SUMIMG+h 

IF(AHS(TERMP) .LT.EPS.AND.ABS(TERMN) .LT. EPS. and. IPaSS.EQ. 1)  GO  To 
1 25 

20  CONTINUE 
25  CONTINUE 

Coo**  COMPUTE  THE  SQUARED  MODULUS  OF  THE  CORRESPONDING  ANALOG 
AGAN=SuPT (2./PI/FS) *STOV 
AM0DSQ=aGAN»»2/  ( 1 . + (OMEG/OMEGS)  » 2 ) 

Cooo*  COMPUTE  THE  SQUARED  MODULUS  AND  PHASE 
SQMODO= (SAVE1o«2+SAVE2o*2)  oGAlNo*2 

SOMOD= (SUMRELoo2*sumIMG**2) ttGAINoo2 
PHASL  = 57. 30ATAN2 ( -SUM I MG . SUMREL ) 

FMOD  = S(JR T (SQMOU) 

Coooo  WRITE  F REOUENCY, SQUARED  MODULUS.MODULUS.  AND  PHASE 
''-RITE  (6,120)  F , AMODSO , SQMODO . SQMOD ♦ FMOD , PHaSE 
120  FORMAT  < 1 H , F 1 5 . 5 ♦ 4E 1 5 . 5 , F 15 . 2 ) 

10  CONTINUE 
GO  TO  11 
1 CONTINUE 
CALL  EXIT 
STOP 
END 


75 
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TIME-SAMP|_tD  F I RS  f-ORDEP  DYNAMIC  SYSTEM  AT  SAMP.  RATE  20  HZ 
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APPENDIX  C 
PROGRAM  TO  EVALUATE 
THE  AUTOSPECTRUM  OF 
SECOND-ORDER  STOCHASTIC  PROCESSES 


Purpose 

This  utility  program  calculates  and  prints  spectral  results  for 
analog  and  digital  implementations  of  a general,  second-order  stochastic 
process.  The  x-process  is  specified  in  terms  of  its  standard  deviation, 
sampling  generation  interval  (for  the  digital  implementation),  corner 
frequency,  and  damping  ratio.  Analog  and  digital  squared  moduli!  are 
calculated  for  the  specified  process  and  a digital  squared  modulus  is 
calculated  for  a Butterworth  process  having  the  specified  corner  frequency. 

A digital  autospectrum  is  calculated  for  the  specified  process  and  an 
analog  autospectrum  is  calculated  for  the  corresponding  Butterworth 
process . 

Description 

The  analog  squared  modulus  is  developed  from  equation  (5.40),  The 
digital  squared  modulus  uses  (5.35)  with  the  coefficients  (a’s  and  b!s) 
developed  from  the  following  generalization  of  equation  (4.6): 

2 

D = go  + 2Cco  + l 
a a 

b„  = (w2  - 2Sw  + 1)/D, 

/ a a 

where  C is  the  damping  ratio.  The  digital  noise  autospectrum,  Tnn,  used 
in  calculating  the  x-process  autospectrum,  Txx>  via  (5.38)  is  obtained 
from  (5.39).  Note  that  the  ratio  of  variances  YXx(0) /Ynn(0)  (4.15a) 

is  a general  result  and  holds  for  any  set  of  the  (a,  b)  coefficients, 
providing,  of  course,  that  the  denominator  in  (4.15a)  is  positive.  The 
squared  modulus  of  the  second-order  digital  Butterworth  process  is  evaluated 
using  equations  (5.36)  and  (5.37).  The  analog  Butterworth  autospectrum 
uses  equations  (1.1)  and  (1.5). 

Input 


The  program  requires:  (a)  a descriptive  alphameric  title  card; 

(b)  the  process  standard  deviation  desired  (in  arbitrary  units);  (c)  the 
sampling  generation  time  interval  (arbitrary  units);  (d)  the  crossover 
frequency  (compatible  units);  (e)  the  damping  ratio;  (f)  the  number 
of  frequencies  at  which  the  spectra  are  to  be  evaluated;  and  (g)  an  array 
of  frequencies  at  which  the  spectra  are  to  be  evaluated  (compatible  units) • 
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Output 


An  example  of  program  output  is  shown  following  the  source  program 
listing.  Output  names  are  defined  below: 


A. GAIN 
D . GAIN 
FREQ 

ANL. SQ.MOD 
DIG. SQ. MOD 
H22 

ANL.SPK 

DIG.SPK 


the  gain  constant  in  the  analog  autospectrum 

the  gain  constant  (noise  autospectum)  in  the  digital  autospectrum 
frequency 

analog  squared  modulus 

Butterworth  digital  squared  modulus 

squared  modulus  of  the  specified  digital  process 

analog  autospectrum 

digital  autospectrum 
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27 
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36 


SJOB  RDH»KP=29»RUN=FREE»LINES=58 

C****  PROGRAM  TO  EVALUATE  THE  AUTOSPECTRUM  OF  SECOND-ORDER 
STOCHASTIC  PROCESSES. 

INPUTS: 

c****  title  = descriptive  alphameric  title 

C****  STUVX  = THE  PROCESS  STANDARD  DEVIATION 

c**«*  t = Sampling  generation  time  interval 

C*«**  FA  = ANALOG  CORNER  FREQUENCY 
damp  = damping  ratio 

C****  MMX  = NUMBER  of  entries  IN  THE  FREQUENCY  ARRAY 

C****  FREQ (M)  = THE  ARRAY  OF  DISCRETE  FREQUENCIES  AT  WHICH  THE 

C*»o*  TRANSFER  FUNCTIONS  ARE  TO  BE  EVALUATED* 

£«««« 

C*»**  THE  ANALOG  SQUARED  MODULUS  IS  ALSO  CALCULATED. 

IMPLICIT  RE  AL*8 ( A-H , O-Z ) 

DIMENSION  FREQ(50) , TITLE (20) 

DATA  T* OP  1/6,28318500/ 

RT2=DSQRT (2.0D0) 

11  CONTINUE 

READ  (5* 100  »END=1 ) TITLE 
100  FORMAT (20A4) 

WRITE  (6,200)  TITLE 
200  FORMAT ( 1H1 ,20A4) 

READ  (5,300)  STDVX , T ,FA , DAMP , MMX 
300  FORMAT (4Fl0.0,I3) 

WRITE  (6,400)  STQVX,T*FA, DAMP, MMX 

40  0 FORMAT  ( 1H0,4X,6HSTDV.X,9X, lHT,8X,2HFA,6X,AhDAMP,2X»8HMX  FREQS/ 

1 1H  » 4F 1 0 • 5 » 7X ♦ I 3 ) 

C#***  READ  IN  FREQUENCIES 

READ  (5,500)  (FREQ(I) ,1=1, MMX) 

500  FORMAT  ( SF  10.0) 

PI=TWOPI/2. 

WA=DTAN(PI#FA*T) 

c****  GENERATE  RHO  — The  RATIO  OF  PROCESS  TO  N0li>E  VARIANCE 
D = WA<f*2  + 2.«DAMP*WA+l . 

A0=WA**2/D 
Al=2.*AO 
A2  = AO 

Bl=2.* (WA**2-1 . ) /D 

B2=  (WA*«2-2.<>UAMP»WA^1  , ) /D 

AA 1=B 1 

AA2=1 .-Bl**2-B2**2 
BB1=1 . +B2 
BB2=-2.«B1<>B2 
CC1=AO*A1*a2* (A1-B1*A0) 

CC2=-2.»a0*a1*B1-2.*a0*A2*B2-2.*a2*B1*  (Al"Bl*AO)  ♦ AO*^*  A1  **2*  A2**2 
RH0=(CC1*BB2-BB1*CC2)/(AA1*BB2-BB1*AA2) 

C****  GENERATE  THE  NOISE  DENSITY 
GAMN=2.»T*STDVX*»2/RH0 
c*»**  generate  The  analog  gain 

AGaIN=2.*Rt2*STDVX**2/PI/FA 
c»«<h>  write  gain  constants 

WRITE  (6,600)  AGAIN, GAMN 

600  FORMAT ( 1H0,4X,6HA. GAIN, 4X,6HD.GAIN/1H  ,2Fl0.5) 

DEVELOP  COEFFICIENTS  for  the  digital  butterworth  MODULUS 

C1=Wa**4 

C2=3»* (Cl+ 1 . ) 
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37 

38 

39 

40 

4] 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 


C3=4.« (Cl-1 . ) 

C4=C1+1. 

AKG1  = TW()PI«T 
ARG2=2 . * AR(j  1 
WRITE.  (6,800) 

C<mmh>  START  FREQUENCY  LOOP 
DO  10  M= 1 , MMX 
F = F RLQ(M) 

COS1=UCOS(F«aRG1 ) 

SlNl=f)bIN(F*AHGl) 

C052=D©0S (F^ARG2) 

SIN2=OSIN(F*AkG2) 

Fft=F/FA 

ANLMDS=l./(  ( 1 .-FR*«2)  »«2  + 4,  *DAMP**2*FR<k>2) 

5PKA  = ANLM[)S<>  AGAIN 

H22=  ( (AO+Al*>CU51*A2<>COS2)<M>2  + ( A 1 »S  I N 1 ♦ A2»S  I M2 ) **2)  / 

1 ( ( 1 .+bl»COSl  +U2*COS?)  *<»2*  ( Bl^S INI ♦H2*S IN2 ) *<>2) 

DIGMDS=H2? 

SPKD  = U I (,MDS*GAMN 

1>IGMDS=C1*  (3.*4.«COS1*COS2)/(C2*C3»COS1*C4k»COS2> 

WRITE  OUTPUT 

WRITE  (6,700)  F , ANLMDS . U I GMDS , H22 * SPK A , SPKO 
700  FORMAT  ( 1 H , 4F  1 2 . 5 , 2D  1 2 . 5 ) 

800  FORMAT  (1H0,8X,4HFREQ,12H  ANL  . SQ . MOD , 12H  DIG.S(1.M0D,9X»3HH22»5x, 

1 7hANL.SPK,5X,7HDIG.SPK) 

10  CONTINUE 
GO  TO  11 
1 CALI.  FAIT 
EMU 
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sample  of  second  order  spectrum  with  corner  freq  = 0.7  hZ  and  damp 
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